In order to produce contractive tension, myosin II couples a chemical cycle of ATP hydrolysis to a mechanical actin–myosin interaction. Even at the simplest level, a minimal chemical cycle involves ATP binding, hydrolysis, and subsequent release of Pi and ADP.^[26,27] Many models for the mechanochemical cycle of single myosin head have been proposed.^[28–31] To study the cooperative interaction between the two myosin heads and the effect of product concentration on SPOC in a simple way, this paper adopts the four-state mechanochemical cycle model for single myosin head, as shown in Fig. 1.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) The four-state mechanochemical model of single myosin head.

In Fig. 1, S is substrate (ATP), D is product (ADP), the associated states of myosin heads to thin filament are shown by squares, and the dissociated states by circles. The number 1 denotes the step of myosin heads binding to the actin filament in the absence of ATP, while its reverse process is marked as number 2. The step 1 may be due to the electrostatic force between myosin heads and actin filaments when the binding sites have been active by Ca²⁺ or ADP. The step 2 partly results from the drag force created by the relative sliding between the thin filament and the thick filament.^[32] The number 3 denotes the step of ATP binding to the myosin heads and the ATP-myosin complex dissociating from thin filament, and the reversal process of this step is usually negligible in the presence of a sufficiently high ATP concentration. The number 4 denotes the step of power stroke, including the states of ATP hydrolysis and Pi release; the number 6 denotes the step of ADP release; the steps 5 and 7 are the reversal process of 4 and 6, respectively.

The mechanochemical cycle models of single myosin head neglect the cooperative behavior between the two heads. But different forms of cooperative behavior have been proposed to interpret the experimental results. One type of the cooperative behavior is positive cooperativity, that is, the binding of one head promotes the attachment of the second head. This positive cooperativity helps to explain the apparently faster attachment rate of cross-bridges during muscle stretch and shortening.^[25] Based on these facts, this paper takes the two heads as an indivisible whole and proposes a mechanochemical model for myosin II dimer, as shown in Fig. 2.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) The ten-state mechanochemical model for myosin II dimer.

In Fig. 2, X_i represents the states of myosin dimer, k_i are the rate constants. Since each head of the myosin II dimer has four states (see Fig. 1), the myosin dimer with two heads shares a total of sixteen states, and the sixteen states can be simplified to ten states. From these ten states, four mechanochemical cycles can be constructed, marked as 1, 2, 3, and 4, respectively, which are the same as the cycles shown in Fig. 1. All state transitions (indicated by dotted lines) belonging to cycle 4 are ignored by this paper because the state X₁₀ is seldom found in the experiments of muscle contraction.^[16] Although there are still various transitions between other states in cycle 4, such as and , these state transitions outside of the mechanochemical cycle are insignificant relative to those state transitions within the cycles.

In this model, the cooperative behavior between the two heads is characterized by the effect of inter-head interaction on the rates in the chemical cycle, and it is expressed quantitatively by the ratio of . Here, the constant k₄ is the binding rate of single myosin head occupied by ATP to thin filament, and k₈ is the association rate of single myosin head while another head has been attached to thin filament. For the mechanochemical cycle model of single myosin head, the value of k₈ is equal to k₄. But in the proposed model for myosin II dimer, the value of k₈ is several times larger than k₄.

According to the mass action law, the set of chemical dynamics equations describing the graph in Fig. 2 can be written as

The adenine nucleotides are transported through skeletal muscle fibers during muscle activation, and the myofilament lattice significantly reduces the diffusion of adenine nucleotides.^[33–35] As a result, from the outside of myofilament to the center of myofilament, the ATP concentration varies from place to place. Due to the reaction–diffusion effect, the concentration of ATP inside myofilament may change with time. Thus, the average concentration of ATP in the local neighborhood of myosin active sites, y, is a time-dependent variable. Since the influence of myofilament lattice on the diffusion process is very complex, it is difficult to establish the reaction–diffusion equation accurately. The focus of this paper is to study the effect of cooperation between two myosin heads on the mechanochemical cycle. Therefore, the detailed analysis of the diffusion process is not conducted, but only a relatively simple mathematical expression is given below.

The diffusion effect of ATP is mainly reflected in the radial direction of muscle fibers. According to Fickʼs law, the diffusion of ATP follows the equation

Assuming that the volume of myofilaments is V, the side area is A, and the average concentration of ATP in myofilament is y, then the amount of ATP diffused into the muscle fiber per unit time is

If the shape of the myofilament is a cylinder with a length of l and a radius of r, then

There is still no experimental data on the concentration gradient of ATP diffusion in the muscle fibers. For the sake of simplicity, we assume that the concentration gradient is equal everywhere inside the muscle fibers, position-independent, but only as a function of the average ATP concentration y. So

In the experiments of SPOC, the concentration of ATP outside the myofilaments is constant. Here, suppose its value is a mM. The ATP concentration in the center of myofilaments is less than the average concentration of ATP due to the reaction–diffusion effect. For simplicity, let us set , where b is the scaling factor between 0 and 1. Thus, due to the diffusion effect, the change in the concentration of ATP per unit time satisfies the following equation:

Considering that the values of and r are in the same order of magnitude, if the diffusion coefficient D is properly chosen, equation (9) can be simplified to

According to Eq. (1), the amount of ATP consumed by myosin heads per unit time should be expressed as

The SPOC is often studied under isotonic and isometric conditions. Yasuda et al. made a deep research into SPOC and got the time-dependent change in the sarcomere length (SL) under the isotonic condition.^[1] To check the proposed mechanochemical model of myosin dimer, the present work establishes the kinetic equation describing the force balance along the myofilament, as shown by

In Eq. (14), the number of power strokes n can be written as

The second term in Eq. (13), , represents the friction force due to the relative sliding of the thin and thick filaments, where η is the fiction coefficient.^[36]

The last two terms in Eq. (13), F_k and F_ex, represent the passive force originated from the elastic element in the myofibril and the external force, respectively. F_ex is a constant at the isotonic condition. According to the literatures,^[37,38] F_k can be expressed as a function of l

To analyze the cooperative effects of the two heads of myosin II dimer, the proportion of myosin II dimer in different states is numerically calculated for given specified rate constants. Table 1 shows the experimentally measured and estimated range of rates for the steps in the mechanochemical cycle and the actual values employed in our model analyses.

Table 1.

Table 1.

Table 1. Model parameters. . Parameter Value Ref. k1 0.1 s−1 [45] k2 0.5 s−1 [39],[40] K3 [40]–[43] k4 22 s−1 [40],[44] k5 10 s−1 [40],[42]–[44] k6 800 s−1 [40],[46] k7 200 s−1 [41],[44] k8 200 s−1 D [34],[35] r [1] a 0.2 mM [1] b 0.1 c 0.5 mM [47] l0 [48],[49] s0 14.34 nm [50] [48],[49] F0 5.5pN [51] η [36] 0.2 [11] β 0.1 [52]

Table 1. Model parameters. .

The constants k₁ and k₂ are related to the level of muscle activation, which depends on the concentration of ADP and Ca²⁺. The value of k₁ is taken as 0.1 s⁻¹, while k₂ changes with various muscle active conditions.^[39,40] The actual value k₂ employed in our model is 0.5 s⁻¹.

The rate k₃ is the first-order rate constant of ATP binding to myosin, and K₃ is the second-order rate constant of ATP binding to myosin. So, , where y is the average concentration of ATP inside the myofilament. The rate constant K₃ is about ,^[40–43] and the concentration of adenine nucleotides is usually reported in units of mM. As a result, the rate K₃ ranges from to .

The constant k₄ is the binding rate of single myosin head occupied by ATP to thin filament, and k₈ is the association rate of single myosin head while another head has been attached to thin filament, as shown in Fig. 2. Some experimental results suggested that the binding of one head promotes the attachment of the second head during muscle stretch and shortening. Therefore, the rate k₈ may be several times of k₄. In this work, the value of k₄ is taken as 22 s^−1[40,44] and the ratio is in the range of 5–10.

The rate constant k₅ is determined by the inverse process of ATP hydrolysis, and the rate of 1 s⁻¹–10 s⁻¹ for k₅ is adopted by this work.^[40,42–44] The constants k₆ and k₇ are related to the reversible steps of ADP release. In the experimental research by Yasuda et al., the concentration of ADP is 4 mM. As the second-order rate constant of ADP binding to myosin (K₂) is about ,^[41,44] the rates of 800 s⁻¹ and 200 s⁻¹ are used for k₆ and k₇, respectively.^[40,46]

The diffusion coefficient of small molecules in solution is around 10⁻⁸ m²/s–10⁻⁹ m²/s, and the radius of myofilament used in experiments is in the ranges of 10⁻⁵ m–10⁻⁶ m. Referring to the related literature,^[1,34,35] the value of ranges from 10 to 100 in our calculation. The concentration of ATP in solution, a, is , while the ratio coefficient b varies between 0.1 and 0.3. The concentration of myosin in skeletal muscle, c, is about 0.15 mM–0.5 mM.^[47]

Substituting the above parameter values into Eqs. (1) and (12), the numerical solution of x_i versus t is given in Fig. 3.

Fig. 3.

Figure Option ViewDownloadNew Window

Fig. 3. (color online) Proportion of myosin II dimer in different states. The values of rate constants are k1 = 0.1 s−1, k2 = 0.5 s−1, , k4 = 22 s−1, k5 = 10 s−1, k6 = 800 s−1, k7 = 200 s−1, k8 = 200 s−1, , a = 0.2 mM, b = 0.1, and c = 0.5 mM. The initial values of xi are x1(0) = 0.105, x2(0) = 0.041, x3(0) = 0.029, x4(0) = 0.061, x5(0) = 0.295, x6(0) = 0.388, x7(0) = 0.021, x8(0) = 0.029, and x9(0) = 0.031; y(0) = 0.039 mM.

The numerical simulation results show that the proportion of myosin II dimers in different states, x_i, periodically changes with time. The muscle contraction force is proportional to the number of power strokes, which is mainly decided by x₅ and x₆ in the given activation level. Thus, the force cyclically changes with time, and this is the basic reason for SPOC. It is also found that all the values of x_i tend to be constants if the rate of k₄ is equal to k₈. This brings about the conclusion that one of the necessary conditions for SPOC is . The significant difference between k₄ and k₈ is the key factor that results in the remarkable nonlinear behavior in the muscle system. This suggests that the inter-head cooperative behavior may play a critical role in muscle contraction.

The simulation results also show that the concentration of the ATP around the local area of active site changes periodically with time, as shown in Fig. 4. In most of the experiments to study SPOC, the concentration of the ATP in solution outside of myofilament was fixed to a certain value. The previous theoretical studies about SPOC, therefore, have never taken the concentration of the ATP as a variable. These theories may ignore the diffusion of the ATP from the solution to myofibril. Considering the diffusion and the consumption of ATP inside myofilament, this paper establishes the reaction–diffusion equation about ATP. The numerical results indicate that the concentration of ATP inside myofilament periodically varies over time. This is a new discovery that needs to be validated.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) The time-dependent change in ATP concentration.

The time-dependent changes in sarcomere length are calculated by Eq. (18). The parameter values employed in Eq. (18) are shown in Table 1. The values of , l₀, and s₀, taken from rabbit psoas, are , , and 14.34 nm, respectively.^[48–50] The force generated by a single myosin head, F₀, is about 5 pN–6 pN.^[51] The external force, F_ex, ranges from 1.83 N/m²–3.27 × 10⁴ N/m² in the experiment by Yasuda et al., corresponding to about 15 pN–25 pN for a single thin filament (the number of which is 1.2 × 10¹⁵/m² in myofibril^[53]). The friction coefficient, η, is taken as in our simulation.^[36] The value of the probability of cross-bridge formation at relaxation, , is taken as 0.2.^[11] Referring to the experimental value of duty ratio,^[52] the value of 0.1 for β is used.

Substituting the above related parameter values into Eq. (18), the numerical simulation of SL versus t at F_ex = 19 pN, 21 pN, and 23.5 pN are given in Fig. 5. The SL–time relationship displays as the saw-tooth waveform, the period of which is about 3 s and does not rely on the external force, while the amplitude of oscillation decreases with the increase of the external force. These results are similar to the experimental observations shown in Fig. 4 of the literature.^[1]

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) The time-dependent sarcomere length during isotonic SPOC under different external loads. The average external loads are (a) 19 pN, (b) 21 pN, and (c) 23.5 pN. The initial values of xi and y are the same as those in Fig. 3.

In the process of muscle contraction, the lattice spacing between the thin and thick filaments varies with the sarcomere length. So the cross-bridge formation probability coupling to the muscle structure is not a constant. If the coupling effect is not taken into account, the value of q will hold as a constant and will never change with the sarcomere length. On this occasion, the numerical simulation shows that the range of SL oscillation will be quite different from the experimental data, or the distortion in wave shape of SPOC may be observed. It implies that the coupling relation between the muscle structure and the probability of cross-bridge formation is an important element for SPOC, but not the intrinsic factor.

In addition, the simulation results show that the term F_k in Eq. (17), the passive force originated from the elastic element, is not essential to SL oscillation, although the value of F_k can influence the amplitude of SL oscillation. This conclusion is consistent with the related experimental observations.

As a further check on our model, we test the dependence of SPOC on the concentration of ADP ([ADP]), that is, the change in SL–time relationship with the parameters of k₁ and k₇. The parameter k₇, the result of multiplying [ADP] and K₂ (the second-order rate constant of ADP binding to myosin), characterizes the influence of [ADP] on the transition processes, , , and . The parameter primarily reflects the effect of [ADP] on the activation level of muscle contraction, but the relationship between and [ADP] is too complicated to give a quantitative analysis about them.

The simulation results show that the wave form, amplitude, and period of SPOC are influenced by [ADP]. Both the amplitude and the period decrease with the increase of k₇, accompanied by the change of oscillation from the saw-tooth waveform to the triangular type. A similar tendency is also observed when k₁ is increased, because the increase in the ADP concentration may also change the activation level of muscle contraction.

The force equation (18) is just suitable to calculate the time-dependent changes in sarcomere length under isotonic condition. Since the intrinsic factor for SPOC is the chemical oscillation induced by the cooperative operation of two myosin heads, it is reasonable to believe that the local SPOC wave may occur in some sarcomeres even when the total length of myofibrils is fixed, and the oscillatory contractions can propagate to the adjacent sarcomere, forming a traveling wave along the myofibril. Therefore, the present cooperative model for myosin dimer can explain the intrinsic reason for SPOC under both isotonic and isometric conditions. But a more precise equation should be established to describe the force balance under the isometric condition.

In our model, the level of muscle activation, the cooperative effect between the two heads of myosin, and the concentrations of ATP and ADP are important conditions for spontaneous oscillation of muscle contraction.

The level of muscle activation is determined by the concentration of Ca²⁺ or ADP in solution and is expressed mainly by the rate constant k₁. When the rate of k₁ is relatively larger or smaller, the muscle is in a contracted or relaxed state. Conversely, only when the value of k₁ is moderate, the sustained oscillations of contractive tension may be generated. Therefore, muscle activation at a moderate level is the basic condition for SPOC.

The cooperative effect between the two heads of myosin dimer is reflected by the ratio of . A large number of simulations show that the steady state solutions of Eq. (1) will be constant and the sarcomere length will not change with time if the double-headed cooperative effect is not taken into account, even though the muscle activation is under a moderate level. On the contrary, if we assume that there is significant cooperativity between the two heads ( ), the periodic solutions can be obtained from the differential equations; thereby the time-dependent change in sarcomere length is also observed. Therefore, in our model, the cooperative mechanism between myosin heads serves as the primary factor for SPOC.

The concentration of ATP and ADP in the solution is the external condition at which the muscle produces spontaneous vibration. The concentration of ATP determines the value of the parameter a, and thus the value of the rate k₃ can be calculated; the concentration of ADP determines the value of k₇ while affecting the value of k₁. The rest of the rate constants k₂, k₄, k₅, k₆, and the value of parameter c, are constant for the given experimental conditions. Some experimental observations show that the concentration ratio of ATP to ADP ([ATP]/[ADP]) has a great influence on SPOC. By the numerical simulations, we find that increasing the concentration of ADP (i.e., increasing the rate of k₇) at a given ATP concentration ([ATP] = a mM) gradually changes the period and the amplitude of SOPC until the vibration disappears. This conclusion is consistent with the observed phenomenon in the experiments.

The reaction–diffusion effect of ATP inside muscle fibers is also a necessary condition for SPOC. Without considering the reaction–diffusion effect, the average value of ATP concentration near the active sites of myosin, y, should be assumed as a constant, then the steady-state solutions of the differential equations are also constant and do not change with time. In this case, the muscle will be in a contracted or relaxed state. However, the reaction–diffusion equation (12) has a periodic solution only if the cooperative effect between myosin heads is fully considered ( ). Therefore, the time-dependent change in the average value of ATP concentration inside the muscle fibers is the result of the double-head cooperative mechanism rather than the cause of SPOC.